Bifurcation, chaos, and stability analysis to the second fractional WBBM model

This manuscript investigates bifurcation, chaos, and stability analysis for a significant model in the research of shallow water waves, known as the second 3D fractional Wazwaz-Benjamin-Bona-Mahony (WBBM) model. The dynamical system for the above-mentioned nonlinear structure is obtained by employing the Galilean transformation to fulfill the research objectives. Subsequent analysis includes planar dynamic systems techniques to investigate bifurcations, chaos, and sensitivities within the model. Our findings reveal diverse features, including quasi-periodic, periodic, and chaotic motion within the governing nonlinear problem. Additionally, diverse soliton structures, like bright solitons, dark solitons, kink waves, and anti-kink waves, are thoroughly explored through visual illustrations. Interestingly, our results highlight the importance of chaos analysis in understanding complex system dynamics, prediction, and stability. Our techniques’ efficiency, conciseness, and effectiveness advance our understanding of this model and suggest broader applications for exploring nonlinear systems. In addition to improving our understanding of shallow water nonlinear dynamics, including waveform features, bifurcation analysis, sensitivity, and stability, this study reveals insights into dynamic properties and wave patterns.


Introduction
Solitons are solitary wave solutions that preserve their structure and amplitude while traveling through a system without dispersing or losing energy due to dispersion or dissipation [1].They are fascinating entities that emerge from the delicate interplay between dispersion and nonlinearity [2], and they have been observed and studied extensively in different disciplines, including fluid dynamics [3], plasma physics [4], condensed matter physics [5], and optics [6].Consequently, the attraction of soliton solutions in the field of mathematical physics is considerable [7,8].Many distinguished nonlinear models provide evidence of the presence of these soliton outcomes like the Jimbo-Miwa structure [9], Phi-4 equation [10], the extended Boussinesq nonlinear framework [11], the Nizhnik-Novikov-Veselov model [12], the ISLW system [13], the 3D fractional WBBM structure [14], Fokas-Lenells system [15], and various others [16][17][18].Nonlinear models can be solved using a variety of effective methods, and soliton outcomes can be obtained.This list includes some techniques, such as the Kudryashov technique [19], the unified process [20], the exp a function approach [21], the modified Sardar subequation algorithm [22], the Hirota bilinear process [23], and various others [24][25][26].
Bifurcation analysis is a robust technique widely used to explore dynamic structures, which has significant implications in various fields [27,28].The bifurcation process, announced by Liu and Li in 2002 [29], is a necessary tool for probing the dynamics of nonlinear models.Particularly efficient at scrutinizing the bifurcation nature and deriving precise soliton outcomes, this method investigates how altering system parameters influences its qualitative behavior [30].By studying bifurcation, investigators gain insights into changes between stable and unstable situations or chaotic dynamics.This manuscript aims to explore novel waveforms and bifurcation analysis of the 2nd fractional WBBM nonlinear structure [31].
The synopsis of the existent investigation is as follows: In Section 2, an explanation of the non-integer order derivative and its basic properties is presented.Section 3 delves into the description and ordinary differential form of the second fractional 3D-WBBM system.Bifurcation analysis is thoroughly described in Section 4, while Section 5 addresses the chaotic activities of the adopted nonlinear structure.The sensitivity examination of the projected nonlinear problem is contained in Section 6, and Section 7 shows the appearance of bright and dark solitons in the model.Section 8 contains the graphical representations conducted for the study.In Section 9, the novelty of the outcomes is addressed, demonstrating original findings attained through the procedures recommended for Eq (7), which have not been previously stated.Finally, Section 10 captures the conclusions drawn from this document.

Conformable derivative and its renowned characteristics
To understand the dynamic processes in function development, the fractional derivative is compared with the numerical derivative, which illustrates their overarching connections.Fractional calculus offers a versatile approach applicable to various domains, such as hydrodynamics, applied mathematics, fluid mechanics, quasi-chaotic dynamical systems, system validation, finance, unpredictable fluid infrastructures, and research methodologies.Additionally, it extends to diverse topics like optical fibers, solid-state biological processes, environmental studies, and theoretical electrical control, among others.A fractional derivative provides a clear explanation for the nonlocal properties of mathematical representations that differ from conventional calculus, which focuses solely on the present state of matter.Over time, various fractional derivatives have been devised to depict significant physical phenomena [32].Riemann-Liouville derivatives modified by Jumarie [33], conformable derivatives of Atangana [34], their beta derivative [35], and derivatives of Caputo [36] provide a better representation than integer-order derivatives.These derivatives find application across diverse domains in contemporary science and engineering.
Consider w : ½0; 1� !R, then the conformable fractional derivative of w(t) with order γ 2 (0, 1] is denoted by D g t wðtÞ and is defined by [37]: Here we include some renowned characteristics of fractional derivatives of specified order γ 2 (0, 1] w. r. to t > 0. If w(t) = w and r(t) = r are any real functions, then (i)

Governing equation
The BBM equation appeared in 1972 as an extension of the KdV model, which was developed to present surface water waves in a homogeneous system.An improvement of the KdV structure, the BBM equation finds applications beyond water surface waves, covering Rossby and drift waves in plasma.The BBM nonlinear system, as outlined in [38], is expressed as follows: and the KdV equation corresponds to the next form: The above Eqs (2) and ( 3) are fundamental instruments for understanding a diverse range of wave properties.They play crucial roles in the study of surface waves in water, acoustic and gravity waves in flexible fluids, hydromagnetic waves in plasmas, nonlinear dispersive processes' long waves, and harmonic crystals' acoustic waves, among other applications.In 2017, Wazwaz introduced a novel equation called the WBBM equation [39], which was derived from a modified three-dimensional BBM equation.This equation is expressed as follows: This article will center on the previously mentioned equation Eq (5) from a fractional perspective, known as the second fractional 3D WBBM equation, as introduced in [31].It is formulated as follows: with real function w(t, x, y, z) of free components t, x, y, and z.D g t ; D g x ; D g y , and D g z implies the non-integer order γ derivatives w.r. to t, x, y, and z, in sequence with 0 < γ � 1, and 0 � t.Applying the next traveling wave relation wðx; y; z; tÞ on Eq (7) with l 1 6 ¼ 0, l 2 6 ¼ 0, l 3 6 ¼ 0, and l 4 6 ¼ 0, one reaches Upon integrating Eq (9) w. r. to B, one arrive at where l 5 implies an integrating constant.For convenience, we set l 5 = 0, then Eq (10) turns into the bellow-mentioned ordinary differential form

Bifurcation analysis
Bifurcation occurs in dynamic systems when small parameter changes cause qualitative changes in the system's behavior [40].It often leads to the appearance of new stable states, periodic orbits, or chaotic behavior.Bifurcation theory helps to understand these sudden changes and predict system behavior in different states [41].This paragraph offers an overview of the bifurcation and phase diagrams of the upcoming planner dynamic framework.By employing this methodology, qualitative analysis of nonlinear models becomes possible.There can be a range of trajectory shapes in this structure, including points, simple closed curves, or similar curves with varying shapes, representing solutions to Eq (7) in a variety of physical structures.
Considering dW dB ¼ P, the planner dynamical framework for Eq (7) can be articulated as follows: By applying the first integral to Eq (12), one reaches the next Hamiltonian function which satisfies the Hamilton canonical equations W 0 ¼ @H @P and P 0 ¼ À @H @W .Here h is an integral constant known as the Hamiltonian constant or energy level.Additionally, it is also called the energy integral or total energy.On the other hand, 1  2 P 2 is the kinetic energy, and a 4 W 4 þ m 2 W 2 is the potential energy of the Hamiltonian system Eq (12).
Consider a real function W(B), which is denoted as a solution to Eq (12) satisfying the physical constraints lim B!À 1 WðBÞ ¼ a 1 and lim B!þ1 WðBÞ ¼ a 2 , with free constants a 1 and a 2 .When a 1 = a 2 , W(B) represents a homoclinic trajectory, yielding w(x, y, z, t) a solitary wave achievement for Eq (7).Moreover, if a 1 6 ¼ a 2 , W(B) represents a heteroclinic orbit.When a 1 > a 2 , w(x, y, z, t) manifests as a kink wave outcome, whereas for a 1 < a 2 , it becomes an antikink wave result for Eq (7).A closed-phase portrait is displayed by Eq (12) in another scenario, which results in a periodic solution for Eq (7).It is worth noting that a phase portrait describes an orbit collection inside a phase plane.
Case 3: α > 0, m > 0 Only one equilibrium point, (0, 0), can be obtained by imposing parameter settings as follows: l 1 = l 3 = 1, and l 2 = l 4 = −1, which is presented in Fig 1c .It is obvious from the figure that (0, 0) denotes a central equilibrium.The paths consist of closed curves that contain a family of periodic orbits (which are red) around (0, 0).As a result, system Eq (12) contains a periodic wave solution.
Case 4: α < 0, m < 0 Only one equilibrium point, (0, 0), can be obtained by imposing parameter settings as follows: It is obvious from the figure that (0, 0) denotes a saddle point.In this case, the trajectories do not encompass closed orbits for the structure Eq (12).

Chaotic behaviors
Chaos refers to a type of behavior exhibited by specific dynamic systems that may appear random but are influenced by deterministic rules [42].Chaotic systems are highly sensitive to initial conditions, meaning even minor alterations in initial values can result in significantly divergent outcomes over time.Although such systems often exhibit intricate and unpredictable behavior, they can also display underlying shapes or structures.Chaotic behavior can occur in a variety of natural and artificial systems, including weather patterns, optics, and fluid dynamics [43,44].This section investigates the chaotic characteristics of the resulting dynamic structure by evaluating the disturbed form.This exploration involves studying 3D and 2D phase views.To begin this examination, here we assume dW dB ¼ P, then the dynamic system is: Here, σ cos(ωt) denotes the disturbed form, where σ represents the intensity and ω denotes the system's frequency.In this paragraph, we investigate how the strength and frequency of the allotment affect the structure described in Eq (15).While keeping the key parameters unchanged (l 3 ), we observe chaotic and quasi-periodic patterns with various frequencies and intensities, as depicted in Figs 2a-2c, 3a-3c and 4a-4c.Fig 2a-2c represents the status of Eq (15) with σ = 0. We depict the system's trajectory position following the allocation's intensity and frequency.As shown in Fig 2, Eq (15) exhibits quasi-periodic behavior in time series plots and 2D and 3D phase diagrams.As shown in Fig 3a-3c, the dynamic structure transitions from a quasi-periodic to a chaotic position with a slight increase in intensity and frequency (σ rises to 0.3 and ω = 2.2).Moreover, the system experiences a chaotic state even for significant disturbances in frequency and intensity (σ rises to 1.4 and ω = 3.9) (see Fig 4a-4c).The multistability of Eq (15) is shown in Fig 5a -5c for , σ = 1.9, and ω = 3.9 with different initial settings.We observe the quasi-periodic nature of Eq (15) for initial value (0.8, 0) and chaotic nature for initial values (0, 0.01) and (0.1, 0.2).Therefore the chaotic behavior of Eq (15) is observed under external disturbances with different initial values.

Sensitivity analysis
Sensitivity analysis assesses how different sources of uncertainty in an input can impact a mathematical model's output [45].It helps to identify which inputs have the most significant impact on output and how variations in inputs are transmitted through the system.Sensitivity analysis helps to understand the strength, reliability, and credibility of model predictions that may guide decision-making in complex systems [46].This segment explores how initial conditions affect the disturbed form described by Eq (15) across various intensities and frequencies while keeping the parameter values unchanged (l ).The resulting time series images for the initial settings (0, 0.01) and (0, 0.02) are depicted in

Soliton waveforms of the adopted nonlinear equation
The purpose of this segment is to scrutinize diverse solitons attained from the stated nonlinear framework using the planner dynamics procedure.
Case 1: α < 0, and m > 0 For h 2 0; À m 2 4a À � , a group of periodic trajectories for the dynamic structure expressed in Eq ( 12) can be obtained.For this case, the Hamiltonian function takes on the following structure: r ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi r ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi with S 1 ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi m 2 þ4ah p a q and S 2 ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi À m a À ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi m 2 þ4ah p a q .Substituting Eq (16) in the initial part of the Hamiltonian framework Eq (12) and then integrating yields: Z Q 0 dR ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi where the integration constant is denoted by B 0 .Hence, we arrive at the following pair of periodic wave results: When we set h ¼ À m 2 4a , the condition S 2 1 ¼ S 2 2 ¼ À m a holds, resulting in the subsequent kink wave outcome if we consider the positive expression and antikink wave outcome if we consider the negative expression: Case 2: α > 0, and m < 0 For h 2 À m 2 4a ; 0 À � , we achieve 2 groups of periodic trajectories of the dynamical system Eq (12).For this condition, the Hamiltonian structure can be expressed in the upcoming form r ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi r ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi whereas S 1 ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi m 2 þ4ah p a q and S 2 ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi À m a À ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi m 2 þ4lh p a q .Substituting Eq (18) in the initial part of the Hamiltonian framework Eq (12) and then integrating yields: Z S 2 W dR ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi and dR ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi where the integration constant is denoted by B 0 .Hence, we arrive at the following pair of periodic wave results: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi If h = 0, we find that S 1 ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi À 2m a p and S 2 = 0.This results in the emergence of the subsequent two solitary wave results: dark soliton if we consider the negative expression and bright soliton if we take into account the positive expression.
then the Hamiltonian function takes on the next structure: r ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi a 2 r ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi whereas S 1 ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi .Substituting Eq (21) in the initial part of the Hamiltonian framework Eq (12) and then integrating yields: Z W 0 dR ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi where the integration constant is denoted by B 0 .Hence, we arrive at the following pair of periodic wave results: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi

Graphical representations
Using appropriate parameter selection, this section explores graphical representations of the acquired results, followed by an explanation of their physical significance.Solution w 4 manifests both bright and dark solitons.To visually present the physical behavior of the specific result w 4 under B 0 ¼ 1; g ¼ 0:5; l 1 ¼ l 2 ¼ 1; l 3 ¼ 1 2 ; l 4 ¼ 2; and h = 0 at y = 1, z = 1, a numerical illustration is depicted in

Stability assessment
To examine the stability of the governing model, the current document will perform the linear stability concept, outlined in a source cited as [47].We will assume the mentioned system's integer order, as specified in equation Eq (4).Subsequently, the perturbed solution will take the below-mentioned form: wðt; x; y; zÞ ¼ a þ x; y; zÞ; ð23Þ whilst a corresponds to a stable-state solution of Eq (4).Combining Eqs ( 23) and ( 4), we achieve Linearizing the current equation as expressed in the form of ν gives, Now, let us presume that the subsequent solution to the present framework Fðt; x; y; zÞ ¼ e iðbxþcyþdzÀ rtÞ ; ð26Þ whilst b, c, d correspond to normalized wave number and ρ corresponds to the frequency of perturbation.Using Eqs ( 26) and (25) yields Solving the present equation Eq (27) for ρ gives rðb; c; dÞ The stability behavior of equation Eq ( 28

Novelty of the outcomes
This paragraph aims to illustrate the uniqueness of our achievements by comparing them with those of previous works.By systematically comparing our results with those in the literature [31,[48][49][50], we highlight the novelty of our results and the uniqueness of our implemented technique.Seadawy and his coauthors obtained periodic and hyperbolic function outcomes to the governing nonlinear problem by employing the simple ansatz technique [48].Mamun et al. extracted trigonometric and hyperbolic function solutions of the mentioned model through the (G 0 /G 2 )-expansion technique [31].Demirbilek derived exponential function solutions for the mentioned model by employing the IBSEF technique [49].Inc and his collaborators presented generalized trigonometric and hyperbolic function solutions to this model through the Sarder-subequation scheme [50].In contrast to these existing procedures and outcomes, our results, denoted as q 1 , q 2 , q 3 , q 4 , and q 5 , diverge significantly, highlighting their novelty.Specifically, we present analyses such as bifurcation, chaos, stability, and sensitivity analysis that were previously ignored.As a result of this study, we provide novel insights into the dynamics and behavior of the model, representing the originality of our contribution to the field.

Conclusion
The bifurcation, chaos, and stability analysis for the second fractional 3D WBBM equation has provided significant insights into shallow water wave dynamics.By applying the Galilean transformation, we achieve a dynamic structure that facilitates a comprehensive analysis of bifurcations.Furthermore, our exploration of various solitary wave findings, encompassing periodic waves, bright solitons, dark solitons, anti-kink waves, and kink waves, provided a nuanced understanding of their properties and existence, as visually showcased through simulations.
As shown in our findings, the integration methods employed throughout our study are effective, concise, and efficient.Our findings not only advance our understanding of nonlinear phenomena in shallow water waves but also suggest the possibility of applying our methodologies to more complex nonlinear systems encountered in contemporary engineering scientific contexts.Overall, this study provides avenues for further exploration and emphasizes the significance of interdisciplinary approaches for addressing complex challenges in modern research.

m 2 þ4ah p a q and S 3 ¼
ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi m a þ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi m 2 þ4ah p a q Fig 7a-7d.As shown in Fig 7a and 7b, a bright waveform is obvious for positive values, while Fig 7c and 7d demonstrate a dark waveform for negative values.Another nature demonstrated in Fig 8a-8d clarifies the physical features of the outcome w 2under B 0 ¼ 1; g ¼ 0:5; l 1 ¼ 2; l 2 ¼ l 3 ¼ À 1; l 4 ¼ 1; and h ¼ 1 4 at y = 1, z = 1.Fig 8a and 8b showcase a kink waveform for positive values, while Fig 8c and 8d expresse an anti-kink soliton for negative values.Outcomes w 1 , w 3 , and w 5 display a periodic waveform.Additionally, to elucidate the physical attributes of the specific solution w 1 under B 0 ¼ 1; g ¼ 0:8; l 1 ¼ 2; l 2 ¼ l 3 ¼ À 1; l 4 ¼ 1; and h ¼ 1 5 at y = 1, z = 1, a graphical representation is offered in Fig 9a and 9b.

Fig 8 .Fig 7 .Fig 9 .
Fig 8. Waveform of solution w 2 : (a,c) cubic plot, (b,d) 2D shape.https://doi.org/10.1371/journal.pone.0307565.g008 ) is plotted in Fig 10a and 10b.It can be observed that equation Eq (28) remains stable propagation for b = 1, c = 1, since for the selected parameters, the figure displays a continuous and smooth function ρ(b, c, d) (refer to Fig 10a).On the other hand, the function ρ(b, c, d) remains unstable propagation for b = 1, d = 1, since for the selected parameters, the denominator bc + 1 approaches zero when c is chosen close to −1 (refer to Fig 10b).